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We present a fully nonlinear study of the development of equilibrium after preheating. Preheating is 
the exponentially rapid transfer of energy from the nearly homogeneous inflaton field to fiuctuations 
of other fields and/or the infiaton itself. This rapid transfer leaves these fields in a highly nonthermal 
state with energy concentrated in infrared modes. We have performed lattice simulations of the 
evolution of interacting scalar fields during and after preheating for a variety of infiationary models. 
We have formulated a set of generic rules that govern the thermalization process in all of these 
models. Notably, we see that once one of the fields is amplified through parametric resonance or 
other mechanisms it rapidly excites other coupled fields to exponentially large occupation numbers. 
These fields quickly acquire nearly thermal spectra in the infrared, which gradually propagates into 
higher momenta. Prior to the formation of total equilibrium, the excited fields group into subsets 
with almost identical characteristics (e.g. group effective temperature). The way fields form into 
these groups and the properties of the groups depend on the couplings between them. We also 
studied the onset of chaos after preheating by calculating the Lyapunov exponent of the scalar 
fields. 
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I. INTRODUCTION 

The theory of inflation has been highly successful in 
explaining many of the initial conditions for the Hot Big 
Bang model as well as providing a mechanism by which 
the seeds of large scale structure were formed. Typi- 
cal models of inflation are based on the slow-roll evolu- 
tion of the homogeneous inflaton scalar field(s) ((>. Infla- 
tion ends when the slow-roll regime is dynamically ter- 
minated and the field(s) begins to oscillate around the 
minimum of its effective potential V{(j)) as in chaotic in- 
flation Q , or "waterfalls" towards the minimum of V as 
in hybrid inflation After inflation the homogeneous 
inflaton fleld(s) decays due to its interactions with other 
fields or its self-interaction. If the inflaton decay into 
other fields were slow as in perturbation theory the cre- 
ated particles would settle into thermal equilibrium as 
the inflaton decayed. However, the decay of the inflaton 
typically occurs via rapid, non-perturbative mechanisms 
collectively known as preheating The character of 
preheating may vary from model to model, e.g. para- 
metric excitation in chaotic inflation ||ic|] and another, 
specific type of preheating in hybrid inflation joj , but its 
distinct feature remains the same: rapid ampliflcation of 
one or more bosonic flelds to exponentially large occupa- 
tion numbers. This ampliflcation is eventually shut down 
by backreaction of the produced fluctuations. The end 
result of the process is a turbulent medium of coupled, 
inhomogeneous, classical waves far from equilibrium 

Despite the development of our understanding of pre- 
heating after inflation, the transition from this stage to 
a hot Friedmann universe in thermal equilibrium has re- 
mained relatively poorly understood. A theory of the 
thermalization of the fields generated from preheating 



is necessary to bridge the gap between infiation and 
the Hot Big Bang. The details of this thermalization 
stage depend on the constituents of the fundamental La- 
grangian C{4>i,Xi, '4>i, A^, h^^,, ...) and their couplings, so 
at first glance it would seem that a description of this 
process would have to be strongly model-dependent. We 
have found, however, that many features of this stage 
seem to hold generically across a wide spectrum of mod- 
els. This fact is understandable because the conditions 
at the end of preheating are generally not qualitatively 
sensitive to the details of inflation. Indeed, at the end of 
preheating and beginning of the turbulent stage (denoted 
by t*), the fields are out of equilibrium. We have exam- 
ined many models and found that at there is not much 
trace of the linear stage of preheating and conditions at 

are not qualitatively sensitive to the details of infia- 
tion. We therefore expect that this second, highly non- 
linear, turbulent stage of preheating may exhibit some 
universal, model-independent features. 

Although a realistic model would include one or more 
Higgs- Yang-Mills sectors, we treat the simpler case of 
interacting scalars. Within this context, however, we 
consider a number of different models including several 
chaotic and hybrid inflation scenarios with a variety of 
couplings between the inflaton and other matter flelds. 

There are many questions about the thermalization 
process that we set out to answer in our work. Could 
the turbulent waves that arise after preheating be de- 
scribed by the theory of (transient) Kolmogorov turbu- 
lence or would they directly approach thermal equilib- 
rium? Could the relaxation time towards equilibrium be 
described by the naive estimate t ~ {naint)~^ , where n 
is a density of scalar particles and ai„t is a cross-section 
of their interaction? If the inflaton cf) were decaying into 
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a field X, what effect would the presence of a decay chan- 
nel (T for the X field have on the thermalization process? 
For that matter, would the presence of a significantly 
alter the preheating of x itself, or even destroy it as 
suggested in [^? How strongly model-dependent is the 
process of thermalization; are there any universal fea- 
tures across different models? Finally there's the ques- 
tion of chaos. It is known that Higgs-Yang-Mills systems 
display chaotic dynamics during thermalization [g| . The 
possibility of chaos in the case of a single, self-interacting 
inflaton was mentioned in passing in , but when we be- 
gan our work it was unclear at what stage of preheating 
chaos might appear, and in what way. 

Because the systems we are studying involve strong, 
nonlinear interactions far from thermal equilibrium it is 
not possible to solve the equations of motion using linear 
analysis in Fourier space. Instead we solve the scalar field 
equations of motion directly in position space using lat- 
tice simulations. These simulations automatically take 
into account all nonlinear effects of scattering and back- 
reaction. Using these numerical results we have been 
able to formulate a set of empirical rules that seem to 
govern thermalization after infiation. These rules quali- 
tatively describe thermalization in a wide variety of mod- 
els. The features of this process are in some cases very 
different from our initial expectations. 

Section gives a brief review of preheating in dif- 
ferent inflationary models. This review should serve to 
motivate our study and place it in the broader context 
of inflationary cosmology. Sections III and IV describe 
the results of our numerical calculations. Section |l| de- 
scribes one simple chaotic inflation model that we chose 
to focus on as a clear illustration of our results, while sec- 
tion IV discusses how the thermalization process occurs 
in a variety of other models. Section ^ describes the on- 
set of chaos during preheating and includes a discussion 
of the measurement and interpretation of the Lyapunov 
exponent in this context. Section VI contains a list of 



empirical rules that we have formulated to describe ther- 
malization after preheating. Section VII discusses these 



results and other aspects of non-equilibrium scalar field 
dynamics. Finally, there is an appendix that describes 
our lattice simulations. 



II. INFLATION AND PREHEATING 

In this section we outline the context where the prob- 
lem of thermalization after inflation arises. In the 
inflationary scenario, the very early universe expands 
(quasi)exponentially due to a vacuum-like equation of 
state. Such an equation of state can arise in a number 
of different ways, most of which are based on a homoge- 
neous condensate of one or more classical scalar fields. 
We will consider two types of inflationary models. The 
first is chaotic inflation ^ with the single scalar field 



potential V {(()). The second is hybrid infiation, which 
involves several scalar fields ||] . The properties of these 
models are widely discussed in the literature. We will 
be dealing only with the decay of the homogenous infia- 
ton condensate into inhomogeneous modes of the same 
or other scalars and the subsequent interactions of these 
inhomogeneous modes as they approach thermal equi- 
librium. Any particles present before or during inflation 
are diluted by the exponential expansion. Thus by the 
end of inflation all energy is contained in the potential 
V{(j)^ ...) of one or more classical, slowly moving, homo- 
geneous inflaton fields. Immediately after inflation the 
background field(s) is moving fast and produces parti- 
cles of the fields coupled to it. These created particles 
are mutually interacting and ultimately must end up in 
thermal equilibrium. However, particles may be created 
so fast that they spend some time in non-equilibrium 
states with very large occupation numbers. 
Consider chaotic inflation with the potential 
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Soon after the end of inflation the homogeneous inflaton 
field (j){t) coherently oscillates around the minimum of 
its potential with an amplitude on the order of a Planck 
mass. The inflaton oscillations decay due to the creation 
of particles interacting with 4>. Let x be another scalar 
field coupling with the inflaton field as \g'^4>'^x^ ■ Parti- 
cles of the X field are produced from the interaction of 
the quantum vacuum state of x with the coherently os- 
cillating classical field The dominant channel for this 
production is the non-perturbative mechanism of para- 
metric excitation. The Xk mode functions exponentially 
increase with time as Xfe — e^*"', where the characteristic 
exponent /i^ is a model-dependent function . The 

copious production of x particles constitutes the first 
stage of preheating after inflation This state can be 
studied with analytical methods developed in ||l0|-p^. 
However, very soon the amplitudes of the inhomogeneous 
modes (i.e. the occupation number nu) of x become so 
large that the back-reaction of created particles must be 
taken into account. The most important back-reaction 
effect will be the rescattering of particles x4' @j 

which is difficult to describe analytically jl^l- Thus, to 
follow the evolution of the interacting scalar fields after 
the first stage of preheating (dominated by parametric 
resonance), one must investigate the full non- linear dy- 
namics of the interacting scalars. 

The Hartree approximation, which is often used for 
problems of nonequilibrium quantum field theory, is in- 
sufficient here for several reasons. It fails when field 
fluctuations have amplitudes comparable with that of 
the background field, which occurs exponentially fast in 
our case. It does not take into account the rescatter- 
ing of particles. Moreover, in the context of preheating 
there are diagrams beyond the Hartree approximation 
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that survive in the N —i- oo hmit and give comparable 
contributions to those included in the Hartree approxi- 
mation 

Fortunately, scalar fields with high occupation num- 
bers can be interpreted as classical waves, and the prob- 
lem can be treated with lattice simulations |^. Such 
simulations provide approximate solutions to nonequi- 
librium quantum field theory problems, and we believe 
they include the leading physical effects. 

Hybrid inflation models involve multiple scalar fields. 
The simplest potential for two-field hybrid infiation is 

+ (2) 

Inflation in this model occurs while the homogeneous 
(f> field slow rolls from large </) towards the bifurcation 
point at (f> = (due to the slight lift of the poten- 

tial in (j) direction). Once (j){t) crosses the bifurcation 
point, the curvature of the a field, = d^V/da^, be- 
comes negative. This negative curvature results in ex- 
ponential growth of a fluctuations. Inflation then ends 
abruptly in a "waterfall" manner. It was recently found 
that there is strong preheating in hybrid inflation, 
but its character is quite different from preheating based 
on parametric resonance. 

One reason to be interested in hybrid inflation is that 
it can be easily implemented in supersymmetric theories. 
In particular, for illustration we will use supersymmetric 
F-term inflation as an example of a hybrid model. 

III. CALCULATIONS IN CHAOTIC INFLATION 

In this section we present the results of our numerical 
lattice simulations of the dynamics of interacting scalars 
after inflation. We discuss in detail one simple model 
that we have chosen to illustrate the general properties 
of thermalization after preheating. The next section will 
discuss thermalization in the context of other models. 



A. Model 

The example we have chosen to focus on is chaotic 
inflation with a quartic inflaton potential. The inflaton 
(j) has a four-legs coupling to another scalar field Xj which 
in turn can couple to one or more other scalars at. The 
potential for this model is 

V=lx<^^ + lgW + lhh'af. (3) 

The equations of motion for the model (^ are given by 

+ 3-</>-^v2</,+ (A02+gV)</' = O (4) 



X + 3-x--V2x+(5'(/.2 + /i2a2)x = (5) 



6, + - ^ + {h\^) a, = 0. (6) 

We also included self-consistently the evolution of the 
scale factor a(<). The model described by these equa- 
tions is a conformal theory, meaning that the expan- 
sion of the universe can be (almost) eliminated from the 
equations of motion by an appropriate choice of vari- 
ables [ pT| . See the appendix for more information on 
the lattice simulations we used to solve these equations, 
including information on the initial conditions and the 
rescaled units we used in the calculations and in the plots 
we show here. 

Preheating in this theory in the absence of the Ui fields 
was described in For ^ X the field x will experi- 
ence parametric amplification, rapidly rising to exponen- 
tially large occupation numbers. In the absence of the 
X field (or for sufficiently small g) 4> will be resonantly 
amplified through its own self-interaction, but this self- 
amplification is much less efficient than the two-fleld in- 
teraction. The results shown here are for A = 9 x 10""'^'' 
(for COBE normalization) and g"^ = 200A. When we 
add a third fleld we use h\ — lOQg^ and when we add a 
fourth field we use /i| — 200.g^. 

B. The Output Variables 

There are a number of ways to illustrate the behavior 
of scalar fields, and different ones are useful for exploring 
different phenomena. The raw data is the value of the 
field f{t,x), or equivalently its Fourier transform fk(t). 
One of the simplest quantities one can extract from these 
values is the variance 

{{fit) ' mf) = j^j d'k\m\'' , (7) 

where the integral does not include the contribution of a 
possible delta function at A: = 0, representing the mean 
value /. 

One of the most interesting variable to calculate is the 
(comoving) number density of particles of the /-field 

nfit) = J^ J d'knkit), (8) 

where Uk is the (comoving) occupation number of parti- 
cles 

n,(t)^^|A.|' + YIM' (9) 
ujk = y^fc2 4- ml^j (10) 
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dp 



For the model 




this effective mass is given by 
hi (ah 



(11) 
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for 0, Xi E^nd (Ti respectively. For the classical waves of / 
that we are dealing with, corresponds to an adiabatic 
invariant of the waves. Formula can be interpreted as 
a particle occupation number in the limit of large ampli- 
tude of the /-field. As we will see below this occupation 
number spectrum contains important information about 
thermalization. Notice that the effective mass of the 
particles depends on the variances of the fields and may 
be significant and time-dependent. The momenta of the 
particles do not necessarily always exceed their masses, 
meaning the interacting scalar waves are not necessarily 
always in the kinetic regime. In particular this means 
that in general we cannot calculate the energies of the 
fields simply as / (Pkukrik because interaction terms be- 
tween fields can be significant. 

From here on we will use n without a subscript to 
denote the total number density for a field, and will use 
the subscript only to specify a particular field, e.g. n^. 
We use ritot to mean the sum of the total number density 
for all fields combined. Occupation number will always 
be written Uk and it should be clear from context which 
field is being referred to. 

In practice it is not very important whether you con- 
sider the spectrum fk and the variance of / or the spec- 
trum rifc and the number density. Both sets of quantities 
qualitatively show the same behavior in the systems we 
are considering. The variance and number density grow 
exponentially during preheating and evolve much more 
slowly during the subsequent stage of turbulence. Most 
of our results are shown in terms of number density n / 
and occupation number because these quantities have 
obvious physical interpretations, at least in certain lim- 
iting cases. We shall occasionally show plots of variance 
for comparison purposes. 

We will follow the evolution of n{t) and nk{t). The 
evolution of the total number density ntot is an indica- 
tion of the physical processes taking place. In the weak 
interaction limit the scattering of classical waves via the 
interaction term ^g^^^X^ can be treated using a pertur- 
bation expansion with respect to g^. The leading four- 
legs diagrams for this interaction corresponds to a two- 
particle collision {(px 4'X)j which conserves ntot- The 
regime where such interactions dominate corresponds to 
"weak turbulence" in the terminology of the theory of 
wave turbulence [|l^. If we see ntot conserved it will 
be an indication that these two-particle collisions consti- 
tute the dominant interaction. Conversely, violation of 
ntot{t) — const will indicate the presence of strong tur- 
bulence, i.e. the importance of many-particle collisions. 



Such higher order interactions may be significant despite 
the smallness of the coupling parameter (and others) 
because of the large occupation numbers nk- Later, when 
these occupation numbers are reduced by rescattering, 
the two-particle collision should become dominant and 
ntot should be conserved. 

For a bosonic field in thermal equilibrium with a tem- 
perature T and a chemical potential /i the spectrum of 
occupation numbers is given by 
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(We use units in which h — 1.) Preheating generates 
large occupation numbers for which equation (|l^) re- 
duces to its classical limit 



T 
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which in turn reduces to nk oc 1/k for k ^ m, jj, and 
nk ~ const, for k <C ni, fjL. We will compare the spectrum 
nk to this form to judge how the fields are thermalizing. 
Here we consider the chemical potential of an interacting 
scalar fields as a free parameter. 

Unless otherwise indicated all of our results are shown 
in comoving coordinates that, in the absence of interac- 
tions, would remain constant as the universe expanded. 
Note also that for most of our discussion we consider 
field spectra only as a function of defined by averag- 
ing over spherical shells in k space. For a Gaussian field 
these spectra contain all the information about the field, 
and even for a non-Gaussian field most useful informa- 
tion is in these averages. This issue is discussed in more 
detail in section ^ 

C. Results 

The key results for this model are shown in Figures [T^ 
which show the evolution of n{t) with time for each 
field and the spectrum nk for each field at a time long 
after the end of preheating. These results are shown for 
runs with one field {(p oiily)j two fields (0 and x), and 
three and four fields (one and two (Ji fields respectively) . 
We will begin by discussing some general features com- 
mon to all of these runs, and then comment on the runs 
individually. 

All of the plots of n{t) show an exponential increase 
during preheating, followed by a gradual decrease that 
asymptotically slows down. See for example Figure 0. 
This exponential increase is a consequence of explosive 
particle production due to parametric resonance. This 
regime is fairly well understood ||l^. After preheating 
the fields enter a turbulent regime, during which n[t) de- 
creases. This initial, fast decrease can be interpreted as a 
consequence of the many-particle interactions discussed 
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above; as nk shifts from low to high momenta the over- 
all number decreases. Realistically, however, the onset 
of weak turbulence should be accompanied by the devel- 
opment of a compensating flow towards infrared modes, 
which we would be unable to see because of our finite 
box size. Thus the continued, slow decrease in n(t) well 
into the weak turbulent regime is presumably a conse- 
quence of the lack of very long wavelength modes in our 
lattice simulations. 




FIG. 1. Number density n for = iA<?!>'' + ^^^(^^X^- The 
plots are, from bottom to top at the right of the figure, n^. 
Tlx, and ntot- The dashed horizontal line is simply for com- 
parison. The end of exponential growth and the beginning 
of turbulence (i.e. the moment t*) occurs around the time 
when ntot reaches its maximum. 

To see why this shift is occurring look at the spectra 
rifc (Figures 16-1^, see also @Jl^). Even long after pre- 
heating Even long after preheating the infrared portions 
of some of these spectra are tilted more sharply than 
would be expected for a thermal distribution (p^ ) . Even 
more importantly, many of them show a cutoff at some 
momentum k, above which the occupation number falls 
off exponentially. Both of these features, the infrared tilt 
and the ultraviolet cutoff, indicate an excess of occupa- 
tion number at low k relative to a thermal distribution. 
This excess occurs because parametric resonance is typ- 
ically most efficient at exciting low momentum modes, 
and becomes completely inefficient above a certain cutoff 
fc* . A clear picture of how the flow to higher momenta 
reduces these features can be seen in Figure ^ which 
shows the evolution of the spectrum Uk for x in the two 
field model. 

Figure || illustrates the initial excitation of modes in 
particular resonance bands, followed by a rapid smooth- 
ing out of the spectrum. The ultraviolet cutoff is ini- 
tially at the momentum A:* where parametric resonance 
shuts down, but over time the cutoff moves to higher k 
as more modes are brought into the quasi-equilibrium 
of the infrared part of the spectrum. Meanwhile the in- 
frared section is gradually flattening as it approaches a 
true thermal distribution. During preheating the exci- 
tation of the infrared modes drives this slope to large, 
negative values. From then on it gradually approaches 
thermal equilibrium (i.e. a slope of —1 to depending 
on the chemical potential and the mass). The relaxation 



time for the equilibrium is significantly shorter than that 
given by formula X/ncjint- This estimate is valid for di- 
lute gases of particles, but in our case the large occupa- 
tion numbers amplify the scattering amplitudes |10||. 
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FIG. 2. Evolution of the spectrum of x in model 
V = ;jA(^'* -I- \g^4>^'X^ ■ Red plots correspond to earlier times 
and blue plots to later ones. For black and white viewing: 
The sparse, lower plots all show early times. In the thick 
bundle of plots higher up the spectrum is rising on the right 
and falling on the left as time progresses. 

Figure ^ shows the evolution of the variances 
((/"/) ) for the two field model. As indicated above it 
shows all the same qualitative features as the evolution 
of n for that model. 
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FIG. 3. Variances for V = \\(t>'^ + \g'^<f>'^x^ ■ The upper 
plot shows {[(f) — ) and the lower plot shows {(x — x)^)- 

We can now go on to point out some differences be- 
tween the models, i.e. between runs with different num- 
bers of fields. The one field model (pure A(/)^) shows the 
basic features discussed above, but the tilt in the spec- 
trum is still very large at the end of the simulation and 
is decreasing very slowly compared to the spectral 
tilt and change in n we see in the two field case. This 
difference occurs because the interactions between cj) and 
X greatly speed up the thermalization of both fields. In 
the one field case (j) can only thcrmalize via its relatively 
weak self-interaction. 

The spectra in the two field run also show a novel fea- 
ture, namely that the spectra for (p and x £^re essentially 
identical, which means among other things 



(15) 



This matching of the two spectra occurs shortly after 
preheating and from then on the two fields evolve iden- 
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tically (except for the remaining homogeneous compo- 
nent of (/)). Aposteriori this result can be understood 
as follows. Looking at the potential X(j)^ + g'^'P'^X^^ the 
second term dominates because of the hierarchy of cou- 
pling strengths = 200A. So the potential V ~ g'^4>^X^ 
is symmetric with respect to the two fields, and therefore 
they act as a single effective field. 

Figures |l^ and |l^ show the effects of adding an ad- 
ditional decay channel for x- The interaction of x ^-iid 
a does not affect the preheating of Xi but does drag a 
exponentially quickly into an excited state. The field a 
is exponentially amplified not by parametric resonance, 
but by its stimulated interactions with the amplified x 
field. Unlike amplification by preheating, this direct de- 
cay nearly conserves particle number, with the result 
that decreases as a grows, and the spectra of <f> and x 
are no longer identical. Instead x ^-nd develop nearly 
identical spectra, 
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and they both thermalize (together) much more rapidly 
than X did in the absence of a. There is a looser rela- 



tionship 



Uy , whose accuracy depends on the 



couplings. The inflaton, meanwhile, thermalizes much 
more slowly; note the low k of the cutoff in the (p spec- 
trum in Figure 18, By contrast, there is no visible cut- 



off in the spectra of x ^nd a and the tilt is relatively 
mild. The most striking property of this chain of interac- 
tion is the grouping of fields; x ^-nd a behave identically 
to each other and differently from (p. This again can 
be understood by the hierarchy of coupling constants, 
/i2 = 100^2 ^ 20, OOOA. The term ft-^^V^ is dominant 
and puts x ^^'^ on an equal footing. 

Varying the coupling h did not change the overall be- 
havior of the system, but it changed the time at which 
a grew. In the limiting case g, a grew with x dur- 
ing preheating and remained indistinguishable from it 
right from the start. (We found this, for example, for 
= 10,00052.) 

When we added a second a field we found that the a 
field most strongly coupled to x would grow very rapidly 
and the more weakly coupled one would then grow rel- 
atively slowly. Note for example that 11^2 in Figure ^ 
grows more slowly than Ua- in Figure 18 despite the fact 



that they have the same coupling to x- In the four field 
case is reduced when the more strongly coupled a 
field grows and this slows the growth of the more weakly 
coupled one. Nonetheless, the addition of another a field 
once again sped up the thermalization of x s-nd the a 
fields. The three fields x, (Ti, and (72 once again have 
identical spectra 
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but in the four field case by the end of the run they 
look indistinguishable from thermal spectra. If there is 



an ultraviolet cutoff for these spectra it is at momenta 
higher than can be seen on the lattice we were using. 
Again, we notice a loose relationship w n^-\-ncri-\-ncr2- 
in this case. 

We close this section with a few words about the 
effective masses of the fields, equation (pT[). All the 
masses are scaled in the comoving frame, i.e. we con- 
sider a^TTig^j, and m is measured in units of momentum 
(see appendix) . Figure ^ shows the evolution of the ef- 
fective masses in the two field model. Note that the 
vertical axis of these plots is in the same comoving units 
as the horizontal (fc) axes of the spectra plots. Since the 
momentum cutoff was of order fc ~ 5 — 10 (see Figure ||) 
the mass of (p was consistently smaller than the typical 
momenta of the field. By contrast started out much 
larger and only gradually decreased. The fluctuations of 
X remained massive through preheating (although with 
a physical mass ^ 1/a) and for quite a while afterwards 
the typical momentum of these fluctuations was k ^ m. 

Field masses 




500 1000 1500 

Effective masses for V 



FIG. 4. , - 

function of time in units of comoving momentum. The lower 
plot is and the upper one is m^. 

Figure |^ shows the evolution of the effective masses for 
the three field model. Once again remains small. Al- 
though mo- grows large briefly it quickly subsides. How- 
ever, m^, with contributions from a and (p, remains rel- 
atively large. Note, however, that the spectrum of x ti^s 
no clear cutoff after a has grown, so it is difficult to say 
whether this mass exceeds a "typical" momentum scale 
or not. 

Field masses 




500 1000 1500 2000 

FIG. 5. Time evolution of the effective masses for the 
model V = jA^"* -f |g^(^^x^ + \h^X^(^^- From bottom to 
top on the right hand side the plots show m^, nicr, and m^. 
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IV. OTHER MODELS OF INFLATION AND 
INTERACTIONS 

The model (^) was chosen to ihustrate our basic re- 
suhs because A(/)^ inflation and preheating is relatively 
simple and well studied. Our main interest, however, 
is in universal features of thermalization. In this section 
we therefore more briefly discuss our results for a variety 
of other models. First we continue with A(/)^ inflation by 
discussing variants on the interaction potential described 
above. Next we discuss thermalization in rn^cf^ models 
of chaotic inflation. Finally we discuss hybrid inflation. 

A. Variations on Chaotic Inflation With a Quartic 
Potential 

We looked at several simple variants of the potential 
. We considered a model with a further decay channel 
for a so that the total potential was 

V = + + \hlx'o' + \hlaW- (18) 

Setting hi = /12 we found that for this four field model 
the evolution of the field fluctuations, spectra, and num- 
ber density were qualitatively similar to those in the four 
field model (||). We found that at late times 

K Tifj ~ n~f < . (19) 

The fields Xi and 7 formed a group with nearly identi- 
cal spectra and evolution and rapid thermalization, while 
(j) remained distinct and thermalized more slowly. Com- 
pare these results to the four field model results in Fig- 
ures ^ and 

We also considered parallel decay channels for 

V = \\^^ + \gl4>\' + \gl4>^^^ + \h\^<y^. (20) 

Setting gi — 32 and = lOOgf we found that at late — 
times 

rici, ~ > K ria- . (21) 

In other words the four fields formed into two groups 
of two, with each group having a characteristic number 
density evolution. 

Finally we looked at adding a self-interaction term for 

X 

V = iA^04 ^ 1^2^2^2 ^ 1 ^^^4 (22) 

with = g^ and found the results were essentially un- 
changed from those of the two field runs with no 
term. The x self-coupling caused the spectra of and 
X to deviate slightly from each other, but their overall 
evolution proceeded very similarly to the case with no x 
self-interaction term. 



B. Chaotic Inflation with a Quadratic Potential 

We also considered chaotic inflation models with an 
rn^c/P inflaton potential. Figures pO[ - p^ show results for 
the model 

y=i™V + i.9W + ^/i\V, (23) 

with TO = IQ-^Mp « 1.22 X 10"GeV" (for COBE), = 
2.5 X IO^toVM^, and h'^ = lOOg^. (See the appendix for 
more details.) We considered separately the case of two 
fields (f) and x and three fields 0, x, and a. This model 
exhibits parametric resonance similar to the resonance in 
quartic inflation | pO| , which results in the rapid growth of 
n seen in these figures. The spectra produced in this way 
are once again tilted towards the infrared. In the two 
field case, (p and x do not have identical spectra as they 
did for quartic inflation. This is because the coupling 
term l/2g^(tP'x^ redshifts more rapidly than the mass 
term 1/2to^</)^, so the latter remains dominant in the 
potential, which is therefore not symmetric between (p 
and X- In the three field case we again see similar spectra 
for X and ct, although they are not as indistinguishable 
as they were in Ac/)'* theory. The basic features of rapid 
growth of n, high occupation of infrared modes, and then 
a flux of number density towards ultraviolet modes and a 
slow decrease in Utot are all present as they were for A0^ 
theory. The shape of the (p spectrum does not appear 
thermal, but it is unclear if this spectrum is compatible 
with Kolmogorov turbulence. 

C. Hybrid Inflation 

Preheating has been studied in many different versions 
of hybrid inflation, mostly only at the early stages when 
the equations for the fluctuations can be linearized. It 
had been thought until recently that preheating was not 
a universal process in hybrid inflation. In our recent 
study 1^, however, we found that there is generally a 
very strong preheating in hybrid models, but its charac- 
ter is quite different from preheating based on parametric 
resonance. We discuss in detail in a separate publica- 
tion 1^ our recent analytical and numerical studies of 
preheating in hybrid inflation models, including a sim- 
ple two-field model (||) as well as more complex SUSY 
F-Term and D-Term models. As with parametric res- 
onance, the result of the instability is the exponential 
growth of long- wavelength modes of the fields. 

In this paper we are mostly interested in preheating in 
the non-inflaton sector and the nonlinear stage after pre- 
heating. In Q we studied the instability in the inflaton 
sector of the hybrid model, i.e. the decay of the homo- 
geneous fields and excitations of their fluctuations. Here 
we take a complementary approach and consider the dy- 
namics of the model with an additional scalar field x 



i) 
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coupled to the fields of the hybrid inflation model. The 
potential is 



V = -|4ES ■ 
4 ' 



,,2|2 



4A|$p (|E| 



(24) 



where A = 2.5 x 10~^ and = 2A. Here $, S and S are 
the complex scalar fields of the inflaton sector and x is 
an additional matter field. Inflation occurs along of the 
<1> direction for ($) > v, when E = E = 0. When the 
magnitude of the slow-rolling field <& reaches the value 
(|<l>c|) = f spontaneous symmetry breaking occurs and 
the E fields become excited. It can be shown that at 
the end of inflation and the start of symmetry breaking 
the complicated potential (p^ ) can be effectively reduced 
to the simple two field potential (|^) (where (j) and a are 
combinations of $ and E, S and — ^A) plus the 
coupling term with /i^x^jEp. 

Figure ^ show the evolution of the six degrees of free- 
dom of the infiaton sector as well as the field x- We see 
that all of the infiaton fields except Im{^) are excited 
very quickly. Later the fields x s-^id /m($) are dragged 
into excited states as well. This dragging corresponds to 
preheating in the non-infiaton sector. The fields x ^md 
Jm(<I>) are excited by their stimulated interactions with 
the rest of the fields. The result of this amplification is 
a turbulent state that evolves towards equilibrium very 
similarly to the chaotic models. Although the details of 
infiation and preheating are very different in hybrid and 
chaotic models, we found that once a matter field has 
been amplified, the thermalization process proceeds in 
the same way. 

Field Variances: F-Term model 




2000 4000 6000 8000 10000 

FIG. 6. Evolution of variances of fields in the model (p4|). 
The two fields that grow at late times, in order of their 
growth, are x smd /m($). 



V. THE ONSET OF CHAOS, LYAPUNOV 
EXPONENTS AND STATISTICS 

Interacting waves of scalar fields constitute a dynam- 
ical system, meaning there is no dissipation and the 
system can be described by a Hamiltonian. Dynami- 
cal chaos is one of the features of wave turbulence. In 



this section we address the question if, how and when 
the onset of chaos takes place after preheating. 

The scalar field fiuctuations produced during preheat- 
ing are generated in squeezed states [|l6|,|o| that are 
characterized by correlations of phases between modes k 
and —k. Because of their large amplitudes we can con- 
sider these fiuctuations to be standing classical waves 
with definite phases. During the linear stage of preheat- 
ing, before interactions between modes becomes signif- 
icant, the evolution of these waves may be or may not 
show chaotic sensitivity to initial conditions. Indeed, 
for wide ranges of coupling parameters parametric res- 
onance has stochastic features 11|, and the issue of 
the numerical stability of parametric resonance has not 
been investigated. When interaction (rescattering) be- 
tween waves becomes important, the waves become de- 
coherent. At this stage the waves have well defined oc- 
cupation numbers but not well defined phases, and the 
random phase approximation can be used to describe the 
system. This transition signals the onset of turbulence, 
following which the system will gradually evolve towards 
thermal equilibrium. 

To investigate the onset of chaos in this system we 
have to follow the time evolution of two initially nearby 
points in the phase space, see e.g. Consider two 

configurations of a scalar field / and /' that are identi- 
cal except for a small difference of the fields at a set of 
points XA- We use f(t,XA), fit,XA) to indicate the un- 
perturbed field amplitude and field velocity at the point 
XA and f'{t, xa), fit, xa) to indicate slightly perturbed 
values at this point. In other words, the field configu- 
rations with f{t, Xa), f{t, Xa) and f'{t, xa), f'{t, xa) are 
initially close points in the field phase space. We then 
independently evolve these two systems (phase space 
points) and observe how the perturbed field values di- 
verge from the unperturbed ones. Chaos can be defined 
as the tendency of such nearby configurations in phase 
space to diverge exponentially over time. This diver- 
gence is parametrized by the Lyapunov exponent for the 
system, defined as 



A.i/o,» 
t ^ Do 



(25) 



where 13 is a distance between two configurations and 
Dq is the initial distance at time 0. Here we define the 
distance D{t) simply as 



Dit)^^y^i\fA~fA\)' + i\f'A~fA 



A 



(26) 



where we define /a = f{t,XA) and the summation is 
taken over all the points where the configurations ini- 
tially differed. 

For illustration we present the calculations for the 
model V — jXc/)^ + \g'^4>'^X^- We did two lattice sim- 
ulations of this model with initial conditions that were 
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identical except that in one of them we multiphed the 
amphtude of x by 1 + 10~^ at 8 evenly spaced points 
on the lattice. Figure |^ shows the Lyapunov exponent 
for both fields 4> and x- Note that the vertical axis is 
\t rather than just A. During the turbulent stage the 
parameter D{t) is artificially saturated to a constant 
because of the limited phase space volume of the sys- 
tem. Fortunately, the most interesting moment around 
ti,, where the chaotic motion begins, is covered by this 
simple approach. Certainly, the field dynamics continue 
to be chaotic in subsequent stages of the turbulence, and 
one can use more sophisticated methods to calculate the 
Lyapunov exponent during these stages H,^. However, 
this issue is less relevant for our study. 



100 200 300 400 500 

FIG. 7. The Lyapunov exponent A for the fields (f} (lower 
curve) and x (upper curve). The vertical axis is \t. 

Both fields show roughly the same rate of growth of A, 
but Ajj. grows much earlier than A^ and therefore reaches 
a higher level. The reason for this is simple. The am- 
plitude of X is initially very small and grows exponen- 
tially, so even in the absence of chaos we would expect 
that during preheating the difference x'i^^ xa) — x(^: ^a) 

must grow exponentially, proportionally to x ~ 

itself. So this exponential growth is not a true indicator 

of chaos. 

To get around this problem and define the onset of 
chaos in the context of preheating more meaningfully we 
introduce a normalized distance function 




E 

A 



J'a ~ Ia 

f'A + fA 



I'a - fA 

Ja + Ia, 



(27) 



that is well regularized even while the field x is being 
amplified exponentially. Figure ^ shows the Lyapunov 
exponent A' = j^og-^^K^ for x- In. this case we see the 
onset of chaos only at the end of preheating. The plot for 
the (j) field is nearly identical. The Lyapunov exponents 
for the fields were A^ k, A'^ ~ 0.2 (in the units of time 
adopted in the simulation). This corresponds to a very 
fast onset of chaos. 




100 200 300 400 500 

FIG. 8. The Lyapunov exponent A' for the fields 
using the normalized distance function A. 



and X 



Thus we see that chaotic turbulence starts abruptly at 
the end of preheating. Initially wave turbulence is strong 
and rescattering does not conserve the total number of 
particles ritot ■ The fastest variation in ntot occurs at the 
same time as the onset of chaos, t* ~ 100 — 200. We 
conjecture that the entropy of the system of interacting 
waves is generated around the moment . As the parti- 
cle occupation number drops, the turbulence will become 
weak and ntot will be conserved. Figure |l| clearly shows 
this evolution of the total number of particles ntot in the 
model. 

We also considered the statistical properties of the in- 
teracting classical waves in the problem. The initial con- 
ditions of our lattice simulations correspond to random 
gaussian noise. In thermal equilibrium, the field velocity 
/ has gaussian statistics, while the field / itself departs 
from that unless it has high occupation numbers. Fig- 
ure ^ shows the probability distribution of the field x 
during the weak turbulence stage after preheating, and 
indeed the distribution is nearly exactly gaussian. Thus, 
at this stage we can treat the superposition of classical 
scalar waves with large occupation numbers and random 
phases as random gaussian fields. 




50 100 150 200 250 

FIG. 9. The probability distribution function for the field 
X after preheating. Dots show a histogram of the field and 
the solid curve shows a best-fit Gaussian. 

During preheating, however, this gaussian distribution 
is altered. A simple measure of the gaussianity of a field 
comes from examining its moments. For a gaussian field 
there is a fixed relationship between the two lowest non- 
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vanishing moments, namely 



(28) 



where S(f> = (f> — {(p) and angle brackets denote ensem- 
ble averages or, equivalently, large spatial averages. We 
measured the ratio of the left and right hand sides of this 
equation for cj) and % and their time derivatives using spa- 
tial averages over the lattice. The results are shown in 
Figures 10 and |Tl|. As expected, the fields are initially 
gaussian, deviate from it during preheating, and rapidly 
return to it afterwards. The plots for the moments of 
the field velocities are similar, although the field veloci- 
ties remain closer to gaussianity. 
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Gaussianity Plot for 




500 1000 1500 2000 

FIG. 10. Deviations from Gaussianity for the field (/> as a 
function of time. The solid, red line shows 3((5(^^)^/(5(j!>*) and 
the dashed, blue line shows 3(50^)^/(50^). 



Gaussianity Plot for x 




500 1000 1500 2000 

FIG. 11. Deviations from Gaussianity for the field x as a 
function of time. The solid, red line shows 3(5x^)^/ i^X^) ^^nd 
the dashed, blue line shows 3(5x^)^/(5x^). 

It is quite important to notice that gaussianity is bro- 
ken around the end of preheating and the beginning of 
the strong turbulence. In particular, it makes invalid the 
use of the Hartree approximation beyond this point. 



VI. RULES OF THERMALIZATION 

This paper is primarily an empirical one. We have 
numerically investigated the processes of preheating and 
thermalization in a variety of models and determined a 



set of rules that seem to hold generically. These rules 
can be formulated as follows: 

1. In many, if not all viable models of inflation there 
exists a mechanism for exponentially amplifying fluctu- 
ations of at least one field x- These mechanisms tend to 
excite long-wavelength excitations, giving rise to a highly 
infrared spectrum. 

The mechanism of parametric resonance in single-field 
models of inflation has been studied for a number of 
years. Contrary to the claims of some authors, this ef- 
fect is quite robust. Adding additional fields (e.g. our a 
fields) or self-couphngs (e.g. x"^) has little or no effect on 
the resonant period. Moreover, in many hybrid models a 
similar effect occurs due to other instabilities. The qual- 
itative features of the fields arising from these processes 
seem to be largely independent of the details of inflation 
or the mechanisms used to produce the flelds. 

2. Exciting one field x is sufficient to rapidly drag all 
other light fields with which x interacts into a similarly 
excited state. 

We have seen this effect when multiple fields are cou- 
pled directly to x ^-^id when chains of fields are coupled 
indirectly to x- il takes is one field being excited 
to rapidly amplify an entire sector of interacting fields. 
These second generation amplified fields will inherit the 
basic features of the x field, i.e. they will have spectra 
with more energy in the infrared than would be expected 
for a thermal distribution. 

3. The excited fields will be grouped into subsets 
with identical characteristics (spectra, occupation num- 
bers, effective temperatures) depending on the coupling 
strengths. 

We have seen this effect in a variety of models. For 
example in the models ^ and (|l^) the x and a fields 
formed such a group. In general, fields that are interact- 
ing in a group such as this will thermalize much more 
quickly than other fields, presumably because they have 
more potential to interact and scatter particles into high 
momentum states. 

4-. Once the fields are amplified, they will approach ther- 
mal equilibrium by scattering energy into higher momen- 
tum modes. 

This process of thermalization involves a slow redis- 
tribution of the particle occupation number as low mo- 
mentum particles are scattered and combined into higher 
momentum modes. The result of this scattering is to de- 
crease the tilt of the infrared portion of the spectrum and 
increase the ultraviolet cutoff of the spectrum. Within 
each field group the evolution proceeds identically for all 
fields, but different groups can thermalize at very differ- 
ent rates. 
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VII. DISCUSSION 

We investigated the dynamics of interacting scalar 
fields during post-inflationary preheating and the devel- 
opment of equilibrium immediately after preheating. We 
used three dimensional lattice simulations to solve the 
non-linear equations of motion of the classical fields. 

There are a number of problems both from the point 
of view of realistic models of early universe preheating 
and from the point of view of non-equilibrium quantum 
field theory that we have not so far addressed. In this 
section we shall discuss some of them. 

Although we considered a series of models of inflation 
and interactions, we mostly restricted ourselves to four- 
legs interactions. (The sole exception was the hybrid 
inflation model, which develops a three-legs interaction 
after symmetry breaking.) This meant we still had a 
residual homogeneous or inhomogeneous inflaton field. 
In realistic models of inflation and preheating we expect 
the complete decay of the inflaton field. (There are rad- 
ical suggestions to use the residuals of the inflaton oscil- 
lations as dark matter or quintessence, but these require 
a great deal of fine tuning.) The problem of residual 
inflaton oscillations can be easily cured by three-legs in- 
teractions. In the scalar sector three-legs interactions 
of the type g'^vcf'x'^ may result in stronger preheating. 
Yukawa couplings h^pc/yip will lead to parametric excita- 
tions of fermions p7| ]. 

There are subtle theoretical issues related to the devel- 
opment of precise thermal equilibrium in quantum and 
classical field theory due to the large number of degrees 
of freedom, see e.g. In our simulations we see the 

flattening of the particle spectra Uk and we describe this 
as an approach to thermal equilibrium, but in light of 
these subtleties we should clarify that we mean approx- 
imate thermal equilibrium. 

Often classical scalar fields in the kinetic regime dis- 
play transient Kolmogorov turbulence, with a cascade 
towards both infrared and ultraviolet modes In 
our systems it appears that the flux towards ultraviolet 
modes is occurring in such a way as to bring the fields 
closer to thermal equilibrium (|lj). Indeed, the slope of 
the spectra Uk at the end of our simulations is close to 
— 1. However, given the size of the box in these simula- 
tions we can little say about the phase space fiux in the 
direction of infrared modes. This question could be ad- 
dressed, for example, with the complementary method 
of chains of interacting oscillators, see This is an 

interesting problem because an out-of-equilibrium bose- 
system of interacting scalars with a conserved number 
of particles can, in principle, develop a bose-condensate. 
It would be interesting to see how the formation of this 
condensate would or would not take place in the con- 
text of preheating in an expanding universe. One highly 
speculative possibility is that a cosmological bose con- 
densate could play the role of a late-time cosmological 



constant. 

The highlights of our study for early universe phe- 
nomenology are the following. The mechanism of pre- 
heating after inflation is rather robust and works for 
many different systems of interacting scalars. There is a 
stage of turbulent classical waves where the initial con- 
ditions for preheating are erased. Initially, before all the 
fields have settled into equilibrium with a uniform tem- 
perature, the reheating temperature may be different in 
different subgroups of fields. The nature of these group- 
ings is determined by the coupling strengths. 
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APPENDIX: THE LATTICE CALCULATIONS 

All of the numerical calculations reported here were 
produced with the program LATTICEEASY, devel- 
oped by Gary Felder and Igor Tkachev. The pro- 
gram and documentation are available on the web at 
http: / /physics. stanford.edu/gfelder/latticeeasy /I . The 
site also includes all the files needed to implement the 
particular models discussed in this paper so anyone can 
easily reproduce our results. In this appendix we merely 
summarize the basics of the calculation; more details can 
be found on the website. All quantities are measured in 
Planck units {Mp w 1.22 x lO^^GeV) and we use / to 
denote a generic scalar field. 

The equations of motion for the fields and the scale 
factor a are solved on a three-dimensional lattice using 
finite differencing for spatial derivatives and a second- 
order staggered leapfrog algorithm for time evolution. 
The evolution equation for a scalar field in an expanding 



universe is 



a 



dV 



(29) 



while the evolution of the scale factor is given by the 
Friedmann equations 



Stt 



'Sp) a 



(30) 



(31) 



11 



where the energy density and pressure of a scalar field 
are given by 



V. 



(32) 



(33) 



In a leapfrog scheme the field values and derivatives are 
known at different times, so it is convenient to combine 
Equations (|3^) and (^ij) to eliminate the field derivatives, 
giving 



(34) 



where the gradient is summed over all fields. 

The initial conditions were set in momentum space 
and then Fourier transformed to give the initial field val- 
ues on the grid. Starting at the end of infiation we gave 
each mode a random phase and a gaussian distributed 
amplitude with rms value 



{\h?) = 



1 



(35) 



expansion, energy was conserved to within half a percent 
over the entire run. 

We also did a number of trials to ensure that our re- 
sults were not sensitive to our time step, box size, or 
number of gridpoints. 

The field equations were simplified by variable redef- 
initions. The redefinitions used and the resulting field 
equations for the chaotic infiation models described in 
the paper are given below. (Details on the hybrid in- 
flation model can be found in [^.) The units for the 
fields, times, and momenta in all the plots in the paper 
are measured in Planck units rescaled as indicated be- 
low. Before these rescalings, time was in physical units 
and distances in comoving coordinates. The momenta k 
are also measured in comoving coordinates and they are 
changed by the rescalings below as 1/x. 



Equations for \<f) 

For the model (^) we redefined the field and spacetime 
variables as 



fpr = —f; Xpr = VA0o^; dtpr = VX(t)0 



dt 
a 



(38) 



where 



dp- 



(36) 



In simulations it's useful to use energy conservation as 
a check of accuracy. Energy conservation in an expand- 
ing universe is described by the equation 



p + 3-(p + p) -0. 
a 



(37) 



In principle one could verify that this equation was be- 
ing satisfied during the run, but in practice p is more 
difficult to evaluate than p. Fortunately there is another 
way to accomplish the same thing. Equatio n (|37| ) can be 
derived from the two Friedmann equations (|3C|) and (|3l|), 
so checking that those two equations are being simulta- 
neously satisfied is equivalent to checking Equation (|37|). 
Since the actual equation for the evolution of the scale 
factor is a combination of these two Friedmann equations 
we were able to check energy conservation by calculating 
the ratio of | to ^p as the program progressed. (We 
verified that checking Equation ( |3l| ) gave the same re- 
sults.) For the A(/)^ runs the theory is nearly conformal, 
so almost the same behavior is obtained with or with- 
out the expansion of the universe (if one uses conformal 
variables). So we duplicated a number of our runs with- 
out expansion and directly checked energy conservation. 
In all cases the results of these two methods of checking 
our accuracy were nearly identical. In every run we did, 
including cases where we did the run with and without 



where 



— .3A2M„ is the value of the inflaton at 



— .•-'T:^^"P 

the end of inflation (i.e. at the start of our simula- 
tions). This value was determined from linear numer- 
ical calculations as the point at which = 0. For 
A = 9 X 10^^"* one unit of program (conformal) time is 
a{\^(j)o)^^tpianck ^ fllO^^^sec and one unit of program 
momentum is VX(j)oEpianck ~ a~^10^^GeV , where 
a is the scale factor. In these variables the evolution 
equations became 



' prVpr 



^pr 



9' 2 



\,pr 



' prXpr H" ^ ^ y^pr 



^ ^i,pr 



0pr = O (39) 



(40) 



— Xpr 



'^i,pr ^ pr'^i,P'r 



a 2 _^ 

A ^P"- a 



^i.pr 



(41) 



a 
a 



87r(/); 



\ fields 

2 ^ r^pr\pr 2 A ^P^'^^^P^ I \ ) 



where primes denote differentiation with respect to ipr 
and angle brackets denote spatial averages over the grid. 
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Equations for m (f) 

For the model (^) we used the foUowing redefinitions 

«3/2 



fp 



-/; Xpr = mx; dtpr — mdt 



(43) 



where in this case (po = .193Mp. For m = 10~^Mp a 
unit of program time corresponded to m~^Tpianck ~ 
10~^^sec and a unit of program momentum corre- 
sponded to a~^mEpianck ~ a~^10^^GeV. The evolution 
equations became 



a V 



pr'rpr 



3a^ 
2 a 



+(bpr + ^4>la-^xlr4'pr = (44) 



Xpr ^ ^ prXpr 



_ _ _ 3a^ 



(45) 



-2v72 3 / a'' 



3a^ 
2"a~ 



+ -^(/)2a ^xL'^».pr = (46) 




[8] 
[9] 
[10 

[11 

[12 
[13 
[14 

[15 

[16 

[17 



[18 
[19 



T. Prokopec and T. Roos, Lattice Study of C lassi- 
cal Inflaton Decay, Phys. Rev. D55, 3768 (1997), ^ep^ 



ph/9610400 



T. Biro, S. Matinyan and B. Miiller, Chaos and Gauge 
Field Theory, World Scientific, Singapore, 1994. 
U. Heinz, C. Hu, S. Leopold, S. Matinyan and B. Miiller, 
Phys. Rev. D56, (1997) 2464. 

L. Kofman, A. Linde, and A. Starobinsky, Towards the 
Theory of R eheating After I nflation, Phys. Rev. D56, 
3258 (1997), |hep-ph/9704452 . 

P. Greene, L. Kofman, A. Linde, and A. Starobinsky, 

Structure of Resonance in Pr eheating After I nflation, 

Phys. Rev. D56, 6175 (1997), |hep-ph/9705347 . 

D. Boyanovsky, H. deVega, R. Holman, D. Lee, A. Singh 

and J. Salgado, Phys. Rev. D54, 7570 (1996). 

G. Felder, L. Kofman, A. Linde, and I . Tkachev, Infla- 

tion After Preheating, JHEP, 8 (2000), tiep-ph/0004024 . 

J. Gacia-Bellido, D. Grigoriev, A. Kusenko and M. 

Shaposhnikov, Phys. Rev. D60, (1999) 123504; hep- 

ph/990244. 

V. Zakharov, V. L'vov, G. Falkovich, Kolmogorov Spec- 
tra of Turbulence, Wave turbulence. Springer- Verlag 
1992. 

D. Polarski and A. Starobinsky, Semiclassicality and De- 
coherence of Cosmologic al Perturbation s, Class. Quant. 
Grav. 13, 377-392, 1996, |gr-qc/9504030|. 



P. Greene and L . Kofman, Phys.Lett. B448 (199 9) 6, 
bep-ph/9807339| ; Phys. Rev. D., 2000, in press; [hep 



ph/0003018 



D. Semikoz and I. Tkachev, Phys.Rev. D55 (1997) 489; 



bep-ph/9507306 



G. Aarts, G . Bonini and C. W etterich, Nucl. Phys. B587 
(2000) 403; [hep-ph/0003262 



[1] A. D. Linde, Particle Physics and Inflationary Cosmology 

(Harwood, Chur, Switzerland, 1990). 
[2] A. Linde, Phys. Lett. B259, 38 (1991); Phys. Rev.D49, 

748 (1994). 

[3] L. Kofman, A. Linde, and A. Starobinsky, Rehe ating 



After Inflat ion, Phys. Rev. Lett. 73, 3195 (1994), |hep- 
th/9405187|. 



[4] S. Yu. Khlebnikov, I. Tkachev, Phys. Rev. Lett. 77, 
(1996) 219 

[5] G. Felder and I. Tkachev, Latticeasy software, 2000, pa- 
per in preparation. 

[6] J. Garcia-Bellido, G. Felder, P. Greene, L. Kofman, 
A. Linde and I. Tkachev, in preparation. 



13 



1 . X lO" 



1 . X lO" 



1. xio" 



NUMBER DENSITY VS. TIME 



Number Density 



500 1000 1500 2000 



FIG. 12. V = l/4:X(t>^. (Note that the vertical scale is 
larger than for the subsequent plots.) 
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FIG. 13. V = l/4,\(jy' + l/2c/20\^ g^/X = 200. The 
upper curve represents n^. 
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FIG. 14. V = 1/4A(^* + \l2g'^(i?^ + l/2/^^x^'^^ 
g^jX = 200, = lOOg^ The highest curve is n^. The 
number density of x diminishes when rio- grows. 
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FIG. 15. V = 1/4A,^^ + l/2ff2<^2^2 ^ i/2hix^ol 
g^/X = 200, hi = 200g^/li = lOOg^. The pattern is sim- 
ilar to the three-field case until the growth of 02- 
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FIG. 16. V = 1/4A(^* 
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FIG. 17. V = 1/4A<;!>'' + l/2g'^ (j>'^x^ , ffV^ = 200. The 
spectra of and x s^re nearly identical. 



Occupation Numiber: t-2400. 




FIG. 18. V = 1/4A(^* + l/2g^(l>^x^ + l/^h'^x^cr^. 
g^ /X = 200, h? — lOOg^ The x ^"^^ spectra are similar, 
but a rises in the infrared). The spectrum of (j) is markedly 
different from the others. 
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FIG. 19. V = 1/4A<^* + 1/232^2^2 ^ i/2ft2^V?, 
g'^/X = 200, hi = 200g^/li = lOOg^ All fields other than 
the inflaton have nearly identical spectra. 
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FIG. 20. V = \l2m^<i? + l/2g^(j)^x 
g^Mp/m? = 2.5 x 10^. The upper curve represents n^. 
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FIG. 21. V = l/2m202 _^ i/2f;202^2 _^ 1/2/12^^2 
g^Ml/rn^ = 2.5 x 10^, = lOOff^. The highest curve i 
Tlx- The field that grows latest is a. 
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FIG. 22. V = l/2m>^ + l/2g^(l>W 
g'^Mp/m'^ — 2.5 x 10^. The upper curve represents the spec- 
trum of X- 
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FIG. 23. V = l/2m>^ + l/2p>^x + l/2/i X , 
g'^Ml/m'^ = 2.5 x 10^ h'^ = lOOg^ The x and a spectra 
are similar, while the spectrum of (j> rises much higher in the 
infrared. 
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